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Abstract. We outline the steps needed in to calibrate the Monte Carlo code in order to per- 
form large scale simulations of real globular clusters. We calibrate the results against A-body 
simulations for N = 2500, 10000 and for the old open cluster M67. The calibration is done by 
choosing appropriate free code parameters. 
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1. Introduction 



o 

Most of the modeling of individual star cluster has been focused on static models based 
on the King model (see Meylan & Heggie 1997). There have been also a small number 
i ' of studies able to follow the dynamical evolution of a system. They were mainly based 
on variants of a Fokker-Planck technique and small A^-body simulations (e.g. Grabhorn 
et al. 1992, Drukier 1993, Murphy et al. 1998, Giersz & Heggie 2003 and Hurley et al. 
2005). 

In this presentation we show the further developments of the Monte Carlo code needed 
l/"") . to properly follow the evolution of real star clusters. The dynamical ingredients of the 
code are basically the same as those described in Giersz (2006 and references therein). 
. This code is based on an original code by Stodolkiewicz (1986), which in turn was based 
7— I ' on the code devised by Henon (1971). The extensions to the code were mainly connected 
with: (i) upgrading the prescription of stellar evolution to the algorithm of Hurley et al. 
(2000), (ii) adding the procedures for the internal evolution of binary stars (Hurley et al. 
. 2002), both using the McScatter interface (Heggie, Portegies Zwart & Hurley 2006), (iii) 
a better treatment of the escape process in the presence of a static tidal field. Then the 
new version of the code was calibrated against the results of A^-body simulations so as 
to allow us to construct dynamical evolutionary models of real star clusters. 



2. The calibration of the Monte Carlo technique 

For Monte Carlo simulations the standard A^-body units are adopted. The unit of 
time, however, is proportional to N/ln(jN), where A^ is a number of stars and 7 is a 
parameter. Additionally, to properly follow the relaxation process two other parameters 
have to be chosen: the range of deflection angles, f3 m in and (3 max = 2f3 m i n and the overall 
time step, r, which has to be a small fraction of the half-mass relaxation time. Properly 
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Figure 1. Evolution of the total mass for TV-body and Monte Carlo models. Left panel for 
different 7, right panel for different iseed - initial random number sequence. Parameters of the 
models are described in the figures. 



chosen above parameters together with additional parameters, which characterize mass 
functions, dynamical interactions of binaries and escape process from the system should 
make it possible to reproduce TV-body simulations and follow the evolution of real star 
clusters. 

For calibrate the Monte Carlo code the TV-body simulations with TV = 2500, 10000 
and 24000 were used. The initial parameters of all TV-body and Monte Carlo runs were 
the same as used by Hurley et al. (2005) for simulations of the old open cluster M67 
(TV = 24000). 

2.1. Models with tidal cutoff 

First, we concentrated on the calibration of Monte Carlo models for which the influence 
of the tidal field of a parent galaxy is characterized by the tidal energy cutoff - all stars 
which have energy larger than E t id c = —GMjr t id are immediately removed from the 
system - G is the constant of gravity, M is the total mass and r t id is the tidal radius. 

As was pointed out by Henon (1975) the value of 7 strongly depends on the mass 
function and distribution of stars in the system. The 7 for equal mass stars is rather well 
known (Giersz & Heggie 1994). For unequal mass case and primordial binaries it is much 
less known. The dependence of results of the Monte Carlo simulations on 7 and initial 
random number sequence is presented in Figure 1. 

The best value of 7 inferred from simulations with TV = 2500 and 10000 is equal to 0.02 
(see Figure 1 left panel). The r and (3 m in are equal to 0.01 and 0.03, respectively. As can 
be seen from Figure 1 (right panel) the spread between models (statistical fluctuations) 
with exactly the same parameters, but with different initial random number sequence 
(iseed) is very substantial. The spread between results with different {3 m in and r is well 
inside the spread connected with different iseed. The statistical fluctuation of the models 
is larger for smaller TV as one can expect. 

2.2. Models with full tidal field 

The process of escape from a cluster for a steady tidal field is extremely complicated. 
Some stars which fulfil the energy criterion (binding energy of a star is greater than 
critical energy Eud s = —1.5(GM/rud)) can be still trapped inside a potential wall. 
Those stars can be scattered back to lower energy before they escape from the system. 
According to the theory presented by Baumgardt (2001) the energy excess of those stars 
is decreasing with the increasing number of stars. So the cluster lifetime does not scale 
linearly with relaxation time as expected from the standard theory. To account for this 
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Figure 2. Left panel: evolution of the total mass for TV-body and Monte Carlo models for 
different a. Right panel: evolution of the number of binaries for TV-body and Monte Carlo 
models for different a. Parameters of the models are described in the figures. 

process in the Monte Carlo code an additional free parameter, a, was introduced. The 
critical energy for escaping stars was approximated by: Eud s — —a[GM/rud), where 
a = 1.5 — a(ln( / yN)/N) 1 ' 4 '. So the effective tidal radius for Monte Carlo simulations is 
r tc// = rud/ct and it is smaller than r t id- This means that for Monte Carlo simulations a 
system is slightly too concentrated comparable to TV-body simulations, but the evolution 
of the total mass is reasonably well reproduced. 

Figure 2 shows the evolution of the total mass and the number of binaries for different 
a for TV = 10000. 

The value of a inferred from the comparison between TV-body and Monte Carlo simula- 
tions for TV = 10000 is equal to about 1.05. The other free parameters for the case of full 
tidal field are the same as for the tidal cutoff case: 7 = 0.02, r = 0.01 and /3 m , n = 0.03. 
Again the spread between models with different /3 TO ,- n and t is well inside the spread con- 
nected with different iseed. The statistical spread also does not substantially influence 
the determination of a. As can be seen on Figure 2. (right panel) the Monte Carlo code 
can reproduce well TV-body simulations not only from respect of the global parameters 
of the system, but also from respect of the properties connected with binary activities. 
Despite the fact that the total number of binaries in the system agrees reasonably well 
with TV-body simulations the total binding energy of the binaries increases too fast for the 
Monte Carlo simulations. This is connected with the fact that the present Monte Carlo 
code does not directly follow the 3- and 4-body interactions as the TV-body code does, 
but uses cross sections. The coalescence of binaries induced by dynamical interactions 
and the exchange interactions are missing in the present Monte Carlo simulations. 

2.3. Model of M67 

The data from TV-body simulations of M67 (Hurley et al. 2005) was used in addition to 
TV = 2500 and 10000 to finally calibrate the Monte Carlo code, namely a. The inferred 
formula is a — 1.5 - 3.0(;n( 7 TV)/TV) 1 / 4 . The comparison of results from TV-body and 
Monte Carlo simulations for M67 confirmed the values of 7, t and (3 m in found for smaller 
TV systems. The results of comparison are summarized in Table 1. Taking into account 
the intrinsic statistical fluctuations of both methods the results presented in the Table 
1 show a reasonably good agreement. At the time of 4 Gyr when the comparison was 
done, both models consists of only a small fraction of the initial number of stars (about 
10%) making the fluctuations even stronger. The Monte Carlo model is slightly too 
concentrated compared to the TV-body one. This can be attributed to the parameter a, 
which leads to smaller effective tidal radius than the tidal radius inferred from TV-body 
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Table 1. Monte Carlo and JV-body results for M67 at 4 Gyr 



I 



JV-body 



MC 



M/Mq 

fb_ 

rt pc 1 
rn pc' 1 

M L /M @ 
M lw /Mq 
rh,Lio pc' 1 



2037 
0.60 
15.2 

3.8 
1488 
1342 

2.7 



1984 
0.59 
15.1 
3.03 
1219 
1205 
2.67 



L - stars with mass above O.5M0 and burning nuclear fuel 
L10 - the same as L but for stars contained within 10 pc 
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Figure 3. Left panel: surface brightness profile for JV-body and Monte Carlo models. Right 
panel: luminosity function for time equal to and 4 Gyr for JV-body and Monte Carlo models. 
Parameters of the models are described in the figures. 
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Figure 4. Left panel: surface density profile for Monte Carlo and JV-body models and for 
observations (Bonatto & Bica 2005). Right panel: luminosity function for the Monte Carlo 
model for time = 4 Gyr and for observations (Montgomery et al. 1993) for main sequence stars 
only, binaries only and both main sequence stars and binaries. Parameters of the models are 
described in the figures. 



simulations. As can be seen from Figure 3 the Monte Carlo code also reproduces well 
the results of JV-body simulations regarding the surface brightness profile and luminosity 
function. 

To finally validate the Monte Carlo model of the old open cluster M67 a brief and 
very preliminary comparison with the observational data (Montgomery et al. 1993 and 
Bonatto & Bica 2005) was performed (Figure 4). Both models: JV-body and Monte Carlo 
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do not reproduce the observations well. They are too centrally concentrated. Also for the 
Monte Carlo model the luminosity function is too shallow for dim stars and too high for 
stars around V = 13 mag. In order to achieve a better agreement with observation the 
initial parameters adopted by Hurley et al. (2005) have to be slightly changed. Definitely, 
more work, simulations and observations are needed. 

3. Conclusions 

It was shown that the Monte Carlo code can be successfully calibrated against small 
.ZV-body simulations. Calibration was done by choosing the free parameters describing 
the relaxation process, such as: coefficient in the Coulomb logarithm 7 = 0.02, minimum 
deflection angle f3 m in, time step r, and coefficient in the formula for the critical energy of 
escaping stars, a = 1.5 — 3.0(Zn(7iV)/A r ) 1 / 4 . The calibrated code successfully reproduced 
the TV-body simulations of the old open cluster M67 (Hurley et al. 2005), which was the 
main objective of the calibration procedure. The code is able to provide as detailed data 
as the observations do. However, it showed also some weaknesses, e.g. some important 
channels of blue stragglers formation are not present (coalescence of binaries due to their 
dynamical interactions) and too crude treatment of the escape process. The work is in 
progress to cure these problems (e.g. a few body direct integrations). It was shown also 
that the Monte Carlo code can be used to model evolution of real star clusters and 
successfully compare results with observations (see Heggie 2008, in this volume). The 
very high speed of the code makes it an ideal tool for getting information about the 
initial parameters of star clusters. It is worth to mention that to complete the model of 
the M67 cluster only about seven minutes are needed! 
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